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Abstract The pace-of-life syndrome (POLS) hypothesis 
predicts that most variation in life history, physiology, 
and behavior among individuals, populations, and 
species falls along a continuum from slow to fast pace 
of life. While there is evidence for climatic gradient- 
mediated POLS patterns among species, this approach 
has rarely been explicitly used to study POLS patterns 
among- and within- populations. In addition, the roles 
of sex in POLS evolution among- or within- populations 
are largely unknown. In this study, we investigated the 
effects of altitudinal gradient and sex on the covariations 
between growth rate and several physiological traits 
closely associated with POLS (blood glucose, baseline- 
and stress-induced glucocorticoids (GCs), hemolysis and 
hemagglutination) in the Asiatic toad Bufo gargarizans. 
Contrary to our expectation, altitudinal gradient had no 
influence on the covariations between growth rate and 
physiological traits, neither at the among- nor within- 
population level, indicating that these trait integrations 
have similar fitness payoffs across hierarchical levels. 
In contrast, we found evidence for sex-specific POLS 
composition: there was a negative covariance structure 
between growth rate and baseline GCs- but only in 
females, and a positive covariance structure between 
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growth rate and baseline GCs- but only in females, and 
a positive covariance structure between growth rate and 
hemagglutination-but only in males. This observation 
indicates that these trait associations differ dramatically 
in advancing fitness for each sex, and supports the idea 
that sex-specific POLS composition could evolve in 
species in which the reproductive roles largely differ 
between the sexes. 


Keywords Bufo gargarizans, constraint, glucocorticoids, 
immunity, metabolism, phenotypic integration, physiological 
pace-of-life syndrome 


1. Introduction 


Life history theory posits that under limiting resources, 
the trade-off between allocation to current versus future 
reproduction or survival can generate a correlative pattern 
between life history traits, resulting in a slow-fast life history 
continuum (Stearns, 1989). Recently, researchers further 
proposed that life-history characteristics and suites of 
physiological (metabolic, immunological and hormonal) and 
behavioral traits have coevolved in response to environmental 
conditions forming a pace-of-life syndrome (POLS) 
(Dammhahn et al., 2018, Montiglio et al., 2018, Careau et al, 
2008, West-Eberhard, 2003, Klingenberg, 2014, Ricklefs and 
Wikelski, 2002, Speakman, 2005, Tomasek et al, 2019). This 


functional integration among multiple phenotypic traits is 
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primarily caused by correlational selection, and such processes 
have important consequences for the evolutionary and 
ecological study of populations (Tieleman et al, 2005, Réale et 
al., 2010). 

Due to global climate change, there is increasing interest 
in the intraspecific local adaptation of populations along 
geographical gradients (Tieleman et al., 2006, Pettersen, 
2020). Populations generally show different life-history 
and physiological traits across latitude/altitude ranges that 
result from distinct selective regimes (Conover et al., 2009). 
Therefore, consistent differences in population trait means 
are anticipated, forming in a POLS across latitude/altitudinal 
gradients. Nevertheless, while there is evidence for latitude/ 
altitude-mediated POLS patterns among species (Tieleman 
et al., 2006; Wiersma et al., 2007, Londono et al., 2015), this 
approach has rarely been explicitly used to study POLS 
patterns among- or within- populations of the same species. 
The few studies looking at differences in within-population 
syndrome structure across different environmental gradients 
revealed large similarities of covariation patterns, however 
these studies were limited to behaviors (Segev et al., 2017; 
Alcalay et al, 2015; Bengston and Dornhaus, 2015; Pruitt et al., 
2008; Debecker and Stoks, 2019; Brans et al., 2018).Considering 
the distinct roles of physiology and behavior in life history 
evolution (Polverino et al., 2018), whether these covariation 
patterns across geographic gradients also hold for physiology 
remains a significant gap in our knowledge. 

Sex-specific optima for reproductive investment and life 
history scheduling could result in sex differences in mean trait 
expression (Wedell et al., 2006). Similarly, as a result of their 
different reproductive roles and the environment, each sex 
may come to differ in the strength of correlation among traits, 
or different traits may covary in males and females. According 
to the source of sex-specific selection on pace of life and trait 
covariances, four conceptual frameworks have been proposed: 
uniform POLS structure, sex-specific POLS with different 
strength trait correlation, distinct POLS structure in each sex, 
no POLS (Hamalainen et al., 2018). For instance, uniform POLS 
would be expected when sexual conflict is either completely 
unresolved- such that high genetic correlations between the 
sexes prevents the evolution of dimorphism, or when sexual 
conflict is absent and therefore the life history strategies of the 
sexes converge. However, classic evidence for these outcomes is 
presently scarce. 

In this study, we examined how metabolic, immunological, 
or hormonal traits covary with one life history trait — growth 
— as well as the role of altitudes and sexes in the among- and 
within populational covariations in the Asiatic toad (Bufo 
gargarizans) (Ricklefs and Wikelski, 2002). The distribution 
range of this species spans an extremely large altitudinal 
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gradient from zero to 4300 m above sea level (Fei, 2009). Their 
life history traits also display significant altitudinal differences 
(Liao and Lu, 2012; Liao et al, 2014; Guo et al, 2016), making 
this species suitable for studies related to climatic gradient- 
mediated or sex-specific POLS. We measured rate of growth, 
baseline blood glucose concentrations (GO), two immunological 
indices (hemagglutination, hemolysis), and glucocorticoids 
(GCs, baseline/stress-induced corticosterone) in a common 
garden setting. As a major blood-borne source of metabolizable 
energy circulating in animal blood, GO may correlate with 
energy-demanding processes such as growth and reproduction 
(Fournier et al, 1992, Tomasek et al., 2019). For innate immunity 
characteristics, we measured hemolysis and agglutination 
since they are more dependent on genetic background than 
adaptive immunity (Tieleman et al., 2005). We measured GC 
concentration because of its functional relationship with 
immunity and energy mobilization (Landys et al, 2006). Many 
studies measuring GC levels have found varia ble associations 
with reproductive traits in fish, amphibians and reptiles 
(Crespi et al., 2013). 

In accordance with the POLS hypothesis and the recent 
POLS-related theoretical framework, we expected a faster 
physiological pace-of-life in B. gargarizans that show higher 
GO and lower GCs and innate immunological indices, and are 
more likely to persist at lower altitudes or in males (Réale et al., 
2010; Hamalainen et al, 2018). We postulated that covariations 
should be formed between growth and physiological traits, 
which are expected to be more pronounced at lower altitudes 
because warmer and more productive environments can 
bring up these associations (Segev et al., 2017; van Noordwijk 
and de Jong, 1986). We also hypothesized that B. gargarizans 
demonstrates sex-specific POLS. However, we were not able 
to predict which sex would display the stronger covariance 
structures due to lack of a priori knowledge on how these 
physiological traits confer to breeding an advantage in this 


species. 
2. Materials and methods 


2.1. Sampling and acclimation procedure A total of 123 
adult males were collected from seven sites in western China 
in summer 2018 and spring 2019. These sites were divided 
into two gradient groups according to their relative altitude 
(group L<~1000m, group H>1000m, Figure S1 and Table 1). 
Sites 1-4 form the first transect of the Min River drainage, 
and sites 5-7 form the second transect along the Dadu River 
drainage. After sampling, all individuals were transferred to 
the animal room in Chengdu Institute of Biology, Chinese 
Academy of Sciences (CIBCAS). A constant temperature 
(20£1)'C and light/dark cycle (12h:12h) was maintained in the 
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animal room. All toads vvere housed individually in a plastic 
container (35.5cmx25cmx15cm, LxWxH) that included a piece 
of wet sponge (5cmx7cm) to preserve humidity and a U shape 
tile (15cmx14cmx7cm) for shelter. The toads were mainly fed 
mealworms (Tenebrio molitor) dusted with calcium powder 
(Exo Terra). To ensure the toads could intake fresh mealworms 
ad libitum, the food plate was cleared and replenished every 
other day. Live crickets (Gryllus rubens) were also supplied once 
a week. We measured the body mass and snout-vent length 
(SVL) once every month. 

After captive acclimation of more than 4 months, each toad 
received intraperitoneal injections of a 50 pL mixture of the 
keyhole limpet hemocyanin (KLH, 2 mg/25 mL, Enzo#ALX- 
202-064-M025) and immunological adjuvants (Sigma#F5881) 
(day 0). Two weeks later, all individuals were injected with 
25 uL of a KLH solution again to further boost the immune 
system (we were not able to acquire those data in the Asiatic 
toads due to the low specification of the designed antibody, 
anti-toad IgY Ab). At day 56, we collected saliva samples for 
a corticosterone assay. At day 63, we sacrificed all the toads 
with 0.2% ms-222 solution and collected plasma by opening 
the abdominal cavity and drawing blood from the opening of 
the right ventricular cavity with a heparin-rinsed syringe. The 
blood sample was centrifuged at 3000 rpm at 4°C for 5 min 
after chilling on ice for <1 h. The upper plasma was carefully 
separated and frozen for further analyses. All experimental 
procedures followed the guidelines of the Animal Care and 
Use Committee of CIBCAS. 


2.2. Baseline and stress-induced corticosterone levels To 
determine the baseline corticosterone content, we collected 
saliva samples. A small, dry cotton ball of known weight 
was inserted in the toad’s mouth for 20 s; all samples were 
collected within 3 minutes in order to minimize the increase in 
corticosterone due to treatment stress (Davis and Maney, 2018). 
To determine stress-induced corticosterone concentrations, 
target toads were placed in cloth bags for 60 min, and the 


Table 1 Sample site and sample size information. 
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second saliva sample was taken thereafter. We examined 
corticosterone content in saliva followed the method proposed 
by (Janin et al, 2012). We weighed the saliva-soaked cotton balls 
and stored them individually at -80°C. To exact corticosterone 
from the saliva samples, we added 1 mL methanol to each 
sample tube for 5 h, transferred the cotton ball and methanol 
to a centrifuge tube with a filter membrane and centrifuged 
at 8000 rpm at 4°C for 5 min. The collected solutions were 
concentrated with a vacuum centrifuge for 4 h at room 
temperature. The dry pellets were dissolved in 60 pL pure 
water for radioimmunoassay. The extracted corticosterone 
samples were measured with Iodine [125I]- Rabbit- Cortisol 
Radioimmunoassay Kit (Beijing North Institute of Biological 
Technology, China). Fetal bovine serum was used as the 
intraplate control. The intra-assay coefficient of variation (CV) 
was<10% and inter-assay CV was<15%. The lower and upper 
effective limits of the corticosterone assay kit were between10 
ng/mL and 500 ng/mL, respectively. 


2.3. Baseline blood glucose content The glucose 
concentration in plasma was determined by a glucose oxidase 
assay kit (Applygen, E1010). Firstly, we diluted the plasma 
samples three folds with PBS buffer, then added 10 pL of the 
standard glucose solution (2000, 1000, 500, 250, 125, 62.5, 31.25, 
15.625 pmol/L) or plasma sample, and 190 pL of the working 
solution to the 96-well plate. After incubation at 37°C for 20 
min, we measured the absorbance value at 550 nm with an 
Enzyme standard instrument (Thermo Scientific Varioskan 
Flash) and calculated the glucose content based on the standard 


curve. 


2.4. Indices of constitutive immunity Complement and 
natural antibody levels were measured using an erythrocyte 
hemolysis- hemagglutination assay (Matson et al., 2005; 
Sparkman and Palacios, 2009). Serial two-fold dilutions of 
50 pL plasma were made with phosphate-buffered saline 
(PBS) in 96-well (8 rowsx12 columns) round (U) bottom assay 
plates. Each well received 50 pL of a rabbit red blood (1%) cell 


Site Gradient Altitude (m) Longitude () Latitude (°) Sample size 
Min River drainage 

1. Yingxiu L 828 103.4532 31.0348 12899 

2. Gengda H 1629 103.3241 31.1058 146139 
3. Wolong H 2053 102.9785 30.863 163142 
4. Dengsheng H 2689 102.3777 30.863 4349 
Dadu River drainage 

5. Shimian IL, 926 102.3777 29.2567 196 

6. Hailuogou H 1689 102.1129 29.6119 133 

7. Kangding H 3239 102.0426 29.8382 58 


[No.1 


suspension. Plates were incubated for 60 min at 20°C and then 
scored. Titers were estimated as the negative log, of the highest 
dilution factor of plasma that showed hemagglutination 
or lysis. Half scores were given for titers that appeared 
intermediate. All the samples were assayed in duplicate with 
positive and negative controls in each plated. 


2.5. Statistical analyses All statistical analyses were performed 
with IBM SPSS Statistics 21.0 software (International Business 
Machines Corporation) and R (Version 4.0.3). We calculated 
the growth rate as the body mass gain per day (g/d) during the 
experiment period. Four subjects died during the experiments 
and were hence excluded from the analysis. To achieve 
normality, blood glucose and corticosterone were transformed 
with Jonson transformation. Erythrocyte hemolysis and 
erythrocyte agglutination were computed with Blom rank- 
based normalization. 

Since our experimental design was not balanced, we 
separated the data into two datasets and analyzed them with 
general linear models. Dataset I included all data on males (data 
from males in transect I and ID), which was used to examine 
the effects of transect and altitude on growth, physiological 
variables and their covariations. Dataset II included the data of 
both sexes from transect IJ, and was used to examine the roles 
of sex and altitude in those phenotypic traits. 

We used a univariate mixed model (UMM) to analyze 
the effects of altitude/transect or altitude/sex on growth 
and physiological variables using the “MCMCglmm” R 
package (Hadfield, 2010). For each UMM, altitude, transect 
(or sex), transect (or sex) and their two-way interaction 
(altitudextransect or altitudexsex) were included as fixed effects, 
SVL or body mass as covariates to control for the effects of age 
(initial body mass but SVL was used for growth rate, since SVL 
and initial mass were correlated and only initial body mass 
was covaried with growth rate). Sampling site was used as a 
random effect. 

To investigate the covariation between growth rate 
and physiological variables at the among- and within- 
population level, we used a bivariate mixed model using the 
“MCMCglmm” R package (Hadfield, 2010). We first examined 
the variance and covariance matrix with covariates in the 
models. We then examined the changes in the variance and 
covariance matrix caused by the inclusion or exclusion of 
fixed factors (altitude, transect/sex) in the bivariate model. 
The response and explanatory variables were scaled before 
analysis (mean centered on 0 and SD reduced to 1) to facilitate 
the interpretation of the results. Sampling site was included 
as a random effect on all traits. We used the following 
priors for the residual (V = diag(2), nu = 0.002) and random 
effect matrices (V = diag(2) x 0.002, nu = 1.002, alpha.mu = 
rep(0.2), alpha.V = diag(25*2, 2, 2)). To compute the posterior 


Ping LI et al. 
Sex Modulate Pace-of-Life Syndrome 


37 


distribution, the model was run over 420 000 iterations, with 
a burn-in of 20 000 and a thinning interval of 100, to obtain 
an effective sample size between 20 001 and 419 901, with an 
autocorrelation level between retained iterations lower than 
0.05. 


3. Results 


Initial body mass was positively associated with growth rate, 
and SVL was negatively associated with hemagglutination 
titer and stress-induced GCs levels (P<0.05, Table S1—S2). 
Moreover, there were no significant differences among 
altitudinal groups in growth rate and physiological variables in 
either dataset (P>0.05, Table S1-S2). Similarly, transect did not 
affect growth rate or physiological indices (P>0.05, Table S1). 
However, growth rate and physiological variables (excluding 
stress-induced GCs) varied between the sexes, with females 
having higher growth rate, hemolysis and hemagglutination 
in females, and males having with higher baseline GCs, higher 
blood glucose content (P<0.05, Figure 1, Ta ble S2). 

At the between-popula tion, we did not find any significant 
covariance between growth rate and physiological variables 
(Table S3-S4). At the within-population level, growth rate 
significantly covaried with baseline GCs in both datasets 
(mean [95% credible interval]: —0.334(—0.610, -0.100] for 
dataset I; -0.568[-0.846, -0.323] for dataset II) (Table S3- 
S4). Moreover, growth rate significantly covaried with 
hemagglutination in both datasets (0.256[0.038, 0.533] for 
dataset I; 0.472/0.230, 0.764] for dataset ID (Table S3-S4). These 
covariances were not modulated by transect or altitudinal 
group (Table S3-S4). Nevertheless, these covariances 
were significantly reduced when sex was included in the 
explained variables of MCMCglmm modeling in dataset 
II (-0.247[-0.426, -0.081] for the combined growth rate and 
baseline GCs, 0.260[0.083, 0.473] for the combined growth rate 
and hemagglutination, Table 2). Further analyses showed 
that the covariance between growth rate and baseline GCs 
was significant in females (male: —0.6(—1.562, 0.921], female: 
-0.225[-0.477, -0.014], Figure 2), and the covariance between 
growth rate and hemagglutination was significant in 
males (male: 0.288[0.0629, 0.565], female: 0.149/—0.047, 0.390], 
Figure 2). 


4. Discussion 


4.1. Effects of altitude and sex on growth and physiological 
traits Our results demonstrate that the life history traits 
and physiological traits closely related with pace-of-life in B. 
gargarizans were similar along the altitudinal gradient. These 
similarities could be largely associated with the mixed effects 
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Figure 1 The effects of sex on snout-vent length (SVL) (A), body mass (B), growth rate (C), blood glucose (D), baseline corticosterone 
(E), stress-induced corticosterone (F), hemolysis (G), hemagglutination (H), of dataset II in Bufo gargarizans. The results are based on 
univariate mixed models (UMM) with altitude, sex, and their interaction as fixed effects, body mass (for growth rate) or SVL (for the other 
dependent variables) as covariates, and sampling site as a random effect. Data were presented as estimated mean values + standard 
error. Stars refer to significant differences between sexes (*P<0.05, **P<0.01). 


of genetic adaptation and phenotypic plasticity (Conover 
and Schultz, 1995). At high altitudes/latitudes, compensatory 
responses in life-history and physiological traits are expected 
to evolve, if the reductions in these phenotypic traits at low 
temperature involves a reduction in fitness, also known 
as counter-gradient variation (Conover and Schultz, 1995; 
Conover et al., 2009). For instance, although tadpoles derived 
from cold environments tend to intrinsically grow faster than 
those derived from warm environments, tadpoles experiencing 
cold conditions at early-life stages grow slower than those 
experiencing warm conditions (Seebacher and Grigaltchik, 
2014). Nevertheless, an earlier study revealed that the bufonids 
exhibit virtually no relation between size at metamorphosis 
and adult size (Werner, 1986). This decoupling suggests that 
maternal effects masking the genetic adaptation in the growth 
of B. gargarizans could emerge after metamorphosis. 

It is noted that in addition to growth rate, the measured 
metabolic, immunological, hormonal indices of B. gargarizans 
differed significantly between the sexes, which could be 
associated with their sex-specific processes of ecology and 
evolution. The high rate of growth in females was traditionally 


considered to be driven by sexual and fecundity selections, 
and substantially contributes to sexual size dimorphism (Chou 
et al., 2016). The relatively low immunocompetence in males 
is likely to be associated with the immunosuppressive roles 
of testosterone content, which synchronously promote the 
development of sexual signals (Folstad and Karter, 1992). In 
addition, the relatively high baseline GCs and blood glucose 
content in males may indicate that the males’ maintenance 
energy expenditure (e.g. respiration, immuno-competency, 
blood circulation, digestion) is higher than females’ (Tomasek 
et al., 2019). 


4.2. Covariation in growth rate and physiological traits 
among and within populations The results of our study 
found strong phenotypic correlations between growth rate, 
and baseline GCs and hemagglutination, supporting POLS 
at the within-population level. However, these correlations 
were not found at the among-population level, indicating 
that the coupling between life-history and physiological traits 
may vary across hierarchical levels (Réale et al., 2010). Since 
high baseline GCs level are suggested to indicate poor health 
condition, the negative correlation between growth rate and 
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Table 2 Estimates of fixed effects and variance components for growth rate and physiological traits in of dataset II, obtained from bivariate 
mixed models. Panel A reports the estimates of fixed effects and variance components only with covariates included in the explanatory 
variables, and panel B reports the variance components after including altitude into the explanatory variables, and panel C reports the 
variance components after including transect into the explanatory variables. Variation in response and explanatory variables were scaled 
before analysis to facilitate the interpretation of the estimates. Given the structure of the data, the residual variance and covariance are 
interpreted as the phenotypic (co)variance within populations. The table gives the mean posterior distribution and its 95% credible interval (CI). 


Growth rate Baseline GCs 
Estimate 95% CI Estimate 95% CI 
Panel A Mar repalzton 0.962 0, 1.421 0.641 0, 1.542 
Var residual 0.997 0.713, 1.313 1.043 0.716, 1.399 
(Cova ies -0.062 -0.545, 0.245 
COVAT residual —0.568 -0.846, —0.323 
Panel B VA alados 21.316 0, 27.827 44.385 0, 48.555 
Var rodas 0.998 0.673, 1.339 1.047 0.700, 1.404 
Covar population -0.763 —4.55, 4.684 
Cova residual -0.568 -0.870, —0.299 
Panel C Mar cute 1.167 0, 4.052 1.188 0, 4.654 
Var residual 0.562 0.377, 0.752 0.813 0.577, 1.094 
Covar population —0.258 —1.447, 1.037 
Covar vesiánal -0.247 -0.426, -0.081 
Growth rate Hemagglutination 
Estimate 95% CI Estimate 95% CI 
Panel A Went riei 0.653 0, 2.354 1.475 0, 4.931 
Val reidil 0.996 0.697, 1.280 0.962 0.652, 1.303 
Covar population 0.083 -0.586, 1092 
Cova pesidual 0.472 0.230, 0.764 
Panel B Mar ropalaton 10.061 0, 42.757 24.355 0, 61.043 
Var eiiiai 0.983 0.691, 1.300 0.93 0.639, 1.244 
Covar population 0.529 -6.136, 8.267 
Cova residual 0.448 0.204, 0.727 
Panel C WEE eren 1.962 0, 6.574 2.597 0, 8.527 
Var residual 0.554 0.363, 0.753 0.863 0.599, 1.160 
Covar population 0.392 -1.352, 3.222 
Covar 0.26 0.083, 0.473 


residual 


baseline GCs could suggest an allocation trade-off between 
growth and survival in the Asiatic toad (Lankford et al, 2001, 
Kim et al., 2011). Moreover, the positive correlation between 
growth rate and hemagglutination level in the Asiatic toad 
is congruent with the prediction of the POLS hypothesis that 
innate immunity is favored by fast pace-of-life (Lochmiller and 
Deerenberg, 2000; Sparkman and Palacios, 2009). 

Although growth rate was integrated with baseline GCs 
and hemagglutination, their covariances did not change 
among or within- altitudinal gradient. This result does not 
support the prediction of POLS that trait integration is 
more likely to occur in warm environments, and these trait 


combinations have similar fitness payoffs (Reznick et al., 2000; 


van Noordwijk and de Jong, 1986). This result is similar to that 
reported in a recent study of Ischnura elegans damselfly larvae, 
which showed that latitude did not alter the covariation 
between life-history and anti-oxidative physiological traits 
(Debecker and Stoks, 2019). On the other hand, a study on the 
water flea Daphnia magan found that the covariation structures 
between life history and stress physiological traits changed 
along an urbanization gradient (Brans et al., 2018). Therefore, 
it seems that the impacts of anthropologic disturbance on 
trait integration already outweigh those imposed by climatic 
gradients. Nevertheless, more investigations in different species 
and manipulation approaches are needed to confirm this 


scenario. 
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Figure 2 The estimate and 95% confidence interval for the covariances between growth rate and baseline GCs and hemagglutination 


for males and females of Dataset II. 


The evolution of sexual dimorphism in POLS has been 
highlighted because of the obvious differences between the 
sexes in life history optima and consequently in the optimal 
expression of life history, behavioral and physiological traits 
involved in POLS (Immonen et al, 2018). In B. gargarizans, we 
detected a negative covariation between growth rate and 
baseline GCs in females and a positive covariation between 
growth rate and hemagglutination in males at the within- 
population level; both of these results indicate a sex-specific 
POLS composition (ie, the specific traits that form a POLS 
differ between the sexes). Considering the obvious sex-biased 
reproductive lifespan in B. gargarizans, these results support 
the predictions that sexual dimorphism in POLS is most likely 
to occur in species whose reproductive roles differ between 
the sexes (Immonen et al., 2018, Hamalainen et al., 2018). The 
enhanced covariation between growth and baseline GCs 
in female toads may help them to invest more resources in 
future reproduction, while the greater covariation between 
growth and hemagglutination in males suggest that they 
are likely to be more sensitive to parasite infections than 
females, which is potentially associated with sex-specific 
hormonal immunosuppression (Desprat et al., 2015). As current 
phenotypic trait covariances may not accurately reflect the 


genetic covariances, further investigations need to address 


the underlying genetic basis of these sex-specific POLS 
composition. 

Overall, this study showed that adult B. gargarizans mainly 
displayed consistent differences in growth rate, metabolic, 
immunological and hormonal traits between sexes, but 
no altitude. Moreover, sex- but not altitude- modified the 
phenotypic integration between growth rate and hormonal 
and immunological traits at the within-population level. 
Our study demonstrates that the prevalence of sex-specific 
POLS and highlights its role in understanding the evolution, 
maintenance and variability of POLS. 
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Figure S1 Sample sites of the study. 
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Table S1 Estimates of fixed effects and variance components for growth rate and physiological traits of dataset I, obtained from univariate 
mixed models. The table gives the mean posterior distribution and its 95% credible interval (CJ). 


Estimate L-95% CI U-95% CI pMCMC 
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Table S2 Estimates of fixed effects and variance components for growth rate and physiological traits of dataset II, obtained from univariate 
mixed models. The table gives the mean posterior distribution and its 95% credible interval (CJ). 


Estimate L-95% CI U-95% CI pMCMC 
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Table S3 Estimates of fixed effects and variance components for growth rate and physiological traits in of dataset I, obtained from bivariate 
mixed models. Panel A reports the estimates of fixed effects and variance components only with covariates included in the explanatory 
variables, and panel B reports the variance components after including altitude into the explanatory variables. Panel C reports the variance 
components after including transect into the explanatory variables. Variation in response and explanatory variables were scaled before 
analysis to facilitate the interpretation of the estimates. Given the structure of the data, the residual variance and covariance are interpreted 
as the phenotypic (co)variance within populations. The table gives the mean posterior distribution and its 95% credible interval (CI). 


Growth rate Blood glucose 
Estimate 95% CI P Estimate 95% CI 


Fixed effects 


"a 


Body mass/SVL -0.211 —0.470, 0.067 0.136 —0.047 —0.291, 0.165 0.633 


Var nt 0.944 0.681, 1.319 0.576 0.419, 0.786 


Covar edat -0.009 -0.202, 0.178 


Var ideat 0.962 0.690, 1.368 0.583 0.390, 0.767 


Covar „esidual -0.01 -0.220, 0.217 
Var residual 0.962 0.683, 1.320 0.58 0.404, 0.788 
Covar sesanal -0.043 -0.272, 0.150 
Fixed effect: Growth rate Baseline GCs 
ixed effects 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL 0.105 -0.163, 0.346 0.462 -0.16 -0.376, 0.101 0.249 


Var punt 0.963 0.669, 1.301 -0.334 -0.610, -0.100 


Covar residual —0.334 -0.610, -0.100 


Variant 0.967 0.673, 1.317 0.939 0.617, 1.248 


Cova sesión -0.335 -0.609, -0.074 


Var aut 0.962 0.678, 1.305 0.91 0.633, 1.202 
Covar residual -0.31 -0.588, -0.078 
I Growth rate Stress—induced GCs 
Fixed effects ——(—————————————————— : 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL —0.197 -0.467, 0.108 0.159 —0.173 —0.431, 0.071 0.2 


Var nal 0.968 0.650, 1.313 0.944 0.651, 1.253 


Covar sexiduar 0.004 -0.266, 0.252 


Var sesidval 0.967 0.691, 1.287 0.958 0.667, 1.307 


Covar residual 0.038 —0.274, 0.300 


Var sesidval 0.958 0.651, 1.289 0.964 0.698, 1.297 


Cova sexiduar 0.009 -0.250, 0.245 


Continued Table S3 


Growth rate Hemolysis 
Estimate 95% CI P Estimate 95% CI P 


Fixed effects 


Body mass/SVL -0.186 0.442, 0.056 0.156 -0.157 -0.418, 0.101 0.279 


Var sesidval 0.981 0.688, 1.349 1.003 0.697, 1.352 


Covar residual 0.155 -0.093, 0.438 


Var sesidval 0.957 0.689, 1.302 1.018 0.704, 1.378 


Covar residual 


Var sión 0.957 0.689, 1.293 0.992 0.687, 1.325 
Covar „esidual 0.122 -0.111, 0.390 
. Growth rate Hemagglutination 
Fixed effects I A 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL -0.181 -0.431, 0.061 0.128 —0.276 —0.576, -0.039 0.046 


Var ideat 0.964 0.681, 1.291 0.881 0.606, 1.193 


Covar residual 0.256 0.038, 0.533 


Var sesidval 0.965 0.682, 1.291 0.866 0.603, 1.205 


Covar sexiduar 0.273 0.046, 0.532 


Var mant 0.944 0.678, 1.262 0.859 0.618, 1.153 


Covar residual 0.245 0.016, 0.497 


Table S4 Estimates of fixed effects and variance components for growth rate and physiological traits in of dataset II, obtained from bivariate 
mixed models. Panel A reports the estimates of fixed effects and variance components only with covariates included in the explanatory 
variables, and panel B reports the variance components after including altitude into the explanatory variables. Panel C reports the variance 
components after including transect into the explanatory variables. Variation in response and explanatory variables were scaled before 
analysis to facilitate the interpretation of the estimates. Given the structure of the data, the residual variance and covariance are interpreted 
as the phenotypic (co)variance within populations. The table gives the mean posterior distribution and its 95% credible interval (CI). 


Growth rate Blood glucose 
Estimate 95% CI Estimate 95% CI 


Fixed effects 


a) 
"a 


Body mass/SVL 0.038 -0.189, 0.293 0.772 -0.206 -0.421, 0.032 0.068 


Var nat 0.995 0.715, 1.351 1.022 0.698, 1.326 


Covar residual 0.091 —0.101, 0.374 


Var ideat 0.99 0.673, 1.300 1.011 0.724, 1.369 


Var residual 0.558 0.382, 0.722 1.016 0.704, 1.335 
Covar sesanal 0.023 -0.131, 0.209 
Fixed effect: Growth rate Baseline GCs 
ixed effects 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL —0.17 —0.392, 0.047 0.146 0.023 —0.204, 0.256 0.849 


Variant 0.997 0.713, 1.313 1.043 0.716, 1.399 


Covar resida -0.568 -0.846, -0.323 


Variant 0.998 0.673, 1.339 1.047 0.700, 1.404 


Covar residual -0.568 -0.870, —0.299 


Var sesidval 0.562 0.377, 0.752 0.813 0.577, 1.094 
Covar residual —0.247 -0.426, -0.081 
. Growth rate Stress—induced GCs 
Fixed effects  -__———  _____—— _ _ _  _ _ ee 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL —0.214 —0.438, 0.040 0.069 -0.436 0.656, -0.195 <0.001* 


Var residual 0.995 0.686, 1.307 0.826 0.594, 1.093 


Covar residual —0.144 -0.366, 0.077 


= 


Var sesidval 0.723, 1.352 0.827 0.583, 1.102 


Covar sexiduar -0.133 -0.390, 0.090 


Var sesidval 0.561 0.387, 0.762 0.838 0.587, 1.034 


Covar pagat -0.101 -0.253, 0.082 


Continued Table S4 


Growth rate Hemolysis 
Estimate 95% CI P Estimate 95% CI P 


Fixed effects 


Body mass/SVL -0.194 0.394, 0.051 0.108 -0.06 -0.305, 0.174 0.664 


Var sesidval 0.981 0.679, 1.286 1.056 0.748, 1.379 
Covar residual 0.232 -0.015, 0.468 
Var sesidval 0.981 0.689, 1.293 1.076 0.755, 1.441 
Covar residual 0.243 0, 0.505 
Var nat 0.557 0.355, 0.730 0.961 0.661, 1.283 
Covar residual 0.025 -0.156, 0.198 
. Growth rate Hemagglutination 
Fixed effects I A 
Estimate 95% CI P Estimate 95% CI P 


Body mass/SVL —0.207 -0.448, 0 0.072 —0.421 0.661, —0.197 <0.001* 


Var sesidval 0.996 0.697, 1.280 0.962 0.652, 1.303 


Covar „eidol 0.472 0.230, 0.764 


Var sesidval 0.983 0.691, 1.300 0.93 0.639, 1.244 


Covar residual 0.448 0.204, 0.727 


Var anal 0.554 0.363, 0.753 0.863 0.599, 1.160 


Covar residual 0.26 0.083, 0.473 


